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ABSTRACT 



The characterization of solar surface differential rotation (SDR) from disk- 
integrated chromospheric measurements has important implications for the study 
of differential rotation and dynamo processes in other stars. Some chromospheric 
lines, such as Ca II K, are very sensitive to the presence of activity on the 
disk and are an ideal choice for investigating SDR in Sun-as-a star observations. 
Past studies indicate that when the activity is low, the determination of Sun's 
differential rotation from integrated-sunlight measurements becomes uncertain. 
However, our study shows that using the proper technique, SDR can be detected 
from these type of measurements even during periods of extended solar minima. 
This paper describes results from the analysis of the temporal variations of Ca II 
K hue profiles observed by the Integrated Sunhght Spectrometer (ISS) during the 
declining phase of Cycle 23 and the rising phase of Cycle 24, and discusses the 
signature of SDR in the power spectra computed from time series of parameters 
derived from these profiles. The described methodology is quite general, and 
could be applied to photometric time series of other Main-Sequence stars for 
detecting differential rotation. 

Subject headings: Sun: activity — Sun: surface magnetism — Sun: rotation — Stars: 
rotation 
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1. Introduction 

Rotation plays a fundamental role in processes of stellar formation and evolution, 
through the modification of hydrostatic balance and the redistribution of chemical elements 
(see Maeder & Eenens 2004 for a review of stellar rotation). The Sun is the only star for 
which surface differential rotation can be studied in great detail, by using spectroscopic 
Doppler- velocity measurements or by tracking the motion of features such as sunspots (Beck 
2000). Furthermore, the internal rotation of the Sun can be probed extensively using the 
techniques developed through helioseismology, that is, the analysis of observed oscillations 
of the solar surface (e.g., Thompson et al. 2003). Although solar surface differential rotation 
(SDR) is still not fully understood, it is thought to be a key ingredient of the global-scale 
solar dynamo that generates the 22-year cycle of magnetic activity (Schrijver & Zwaan 
2000). 

Stellar differential rotation cannot be observed directly since the surface of stars 
usually cannot be imaged with sufficient resolution. However, one may expect active stars 
to be rotating differentially assuming that Sun-like dynamos also work on other stars. 
Several photometric and spectroscopic techniques have been developed to investigate the 
rotation of distant stars. Photometric methods are designed to detect the signature of 
spots in stellar hght-curves. If a star is rotating differentially, large starspots located at 
different latitudes may produce modulations in the light-curve signal that can lead to the 
detection of differential rotation (e.g., Frohlich et al. 2012; Frasca et al. 2011). Flux 
measurements in the calcium H and K spectral lines have also been successfully employed 
to monitor chromospheric stellar activity and detect rotation (e.g., Donahue et al. 1996). 
Spectroscopic techniques measure the change in wavelength position and line profile of 
rotationally-broadened absorption lines produced by the spot-observer relative motion 
due to stellar rotation (e.g., Barnes et al. 2005; Ammler-von Eiff & Reiners 2012). The 
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deployment of asteroseismic space projects, such as the NASA Kepler mission (Koch et al. 
2010), make the detection of stellar differential rotation even more accessible. The Kepler 
satellite provides high-quality optical photometric light-curves of unprecedented precision 
and duration (Borucki et al. 2010; Koch et al. 2010) that routinely allow detailed studies 
of stellar magnetic activity on late-type stars. Rotational modulation due to starspots is 
commonly seen in the Kepler light-curves of these stars, allowing detailed study of the 
surface distribution of their photospheric magnetic activity (e.g., Frohlich et al. 2012; 
Prasca et al. 2011; Savanov & Dmitrienko 2012; Walkowicz & Basri 2011). The multi-year 
duration light-curves provided by the Kepler mission will make it possible to investigate 
in detail how activity phenomena ~ such as starspot evolution, differential rotation, and 
activity cycle - operate on stars with a wide range of masses and convection zone depths. 

Assessing the properties of SDR from disk-integrated measurements has a very 
important diagnostic value for the interpretation of stellar observations. Since the Sun can 
be observed in both disk-integrated and disk-resolved modes, solar observations can serve 
as a test-bed for different techniques that can be applied to other stars in search for their 
differential rotation. Here we present a very effective method for extracting the rotational 
components from time series of photometric measurements. We have applied this method 
to solar disk-integrated chromospheric and photospheric time series, allowing us to compare 
the detected rotational periods to the well known profile of SDR. 

Different solar activity indices and methods have been used in previous studies in 
order to unambiguously extract the signature of SDR from Sun-as-a-star data. Hasler et al. 
(2002) used UARS SOLSTICE data and a combination of wavelet and Fourier analyses, 
but beside the basic period of ~27 days their results were inconclusive. Hempelmann & 
Donahue (1997) applied wavelet analysis to disk-integrated Ca 11 K line core emission time 
series covering the duration of Cycle 21 and 22. Their results match the SDR profile only 
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qualitatively. In a later study, Hempelmann (2002) used a similar time-frequency analysis 
on SOLSTICE UARS data and found a clear signature of SDR from the end of Cycle 22 
until the beginning of Cycle 23. However, SDR was not visible in the foregoing part of 
Cycle 22. Attempts have been made using white-light measurements, but the results are 
negative (Hempelmann 2002). Donahue & Keil (1995) reported a successful detection of 
SDR from the analysis of chromospheric full-disk Ca II K observations made almost daily, 
since 1977, by the National Solar Observatory at Sac Peak K-line monitoring program (later 
referred to as SPO data, Keil & Worden 1984). Donahue & Keil divided the time series in 
200-day intervals and computed the power spectrum of the individual segments using the 
Lomb-Scargle periodogram. Figure 3 of Donahue & Keil (1995) shows the observed solar 
period during Cycles 21 and 22. However, this result was criticized by Schrijver (1996) since 
the derived sidereal rotation period of 28.5 days at the beginning of Cycle 22 is too large 
and the observed decline in rotation period during that cycle is too steep. An additional 
limitation in the analysis is that the Lomb-Scargle periodogram typically produces a large 
number of peaks spread over a broad spectral band. In their investigation, however, only 
the periodicity corresponding to the location of the peak with maximum power within 
a small spectral window centered at around 27 days was considered. Peaks outside that 
frequency range were ignored. 

Ideally, the power spectrum of time series containing the possible signature of SDR 
should have a minimum number of peaks to facilitate the interpretation of the different 
components of the signal. This is particularly important when prior knowledge of the actual 
rotation period is not available, e.g., in the case of stellar differential rotation. Our choice 
of spectral estimator was made to specifically address this point. The detection of stellar 
rotation hinges on the presence of active features (e.g., starspots), and properties such as 
lifetime and emergence rate of these features may play a critical role in the interpretation of 
the observed signal. In the next section we describe the main properties of the time series 
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investigated in this study. In Section 3 we delineate the main steps of our data analysis, 
with particular emphasis on the choice of spectral estimator and computation of statistical 
significance of the derived estimates. The results of the data analysis are discussed in 
Section 4, which also includes the description of a simple model to estimate the effects of 
magnetic field diffusion and the rate of sunspot emergence on the detection of solar rotation 
in Sun-as-a-star observations. 

2. The ISS Ca II K parameter time series 

In this study we use observations from the Integrated Sunlight Spectrometer (ISS), one 
of three instruments comprising the Synoptic Optical Long-term Investigations of the Sun 
(SOLIS, Balasubramaniam & Pevtsov 2011 and reference therein). The ISS, in operation 
since late 2006, takes high spectral resolution (R = 300,000) daily observations of the 
Sun-as-a-star in nine different wavelength bands. For this study we select the observations 
taken in the Ca II K spectral line centered at 393.37 nm, the longest continuous data in the 
ISS data set. A least-square fit to known wavelength positions of five photospheric spectral 
lines situated in the red and blue wings of the Ca II K line is performed to determine the 
linear dispersion. The error in determination of dispersion is extremely small to the point 
that one can see slight 0.02% annual variations that are attributed to seasonal changes in 
the refraction index of air (Bertello et al. 2011). In late 2011, the ISS was upgraded to 
include optical encoders on the spectral grating and mechanisms moving the CCD camera 
(for fiatfielding purposes). This upgrade further improved the ISS stabihty. 

The spectral band for the Ca II K line covers about 0.05 nm, and thus, does not include 
the continuum near the Ca K line. To overcome this limitation, the observed line profiles 
are normalized using intensities at two narrow bands situated in the blue (393.147-393.153 
nm) and red (393.480-393.500 nm) wing of the Ca K line. Mean intensities in these two 
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bands are scaled to match intensities in the "reference" spectral line profile taken by the 
NSO Fourier Transform Spectrometer (FTS, Wallace et al. 2007). This normalization 
also helps to remove any residual linear gradients in intensity in the spectral direction. 
The normalization bands are situated far in the wings of the K line, and thus, are not 
likely affected by solar cycle variation. As a test, we examined intensities in the cores of 
photospheric lines situated near the bands used for scaling the intensities, and did not find 
any systematic trends indicating a solar cycle dependency. For a more detailed description 
of the ISS data reduction we refer the reader to Bertello et al. (2011). 

From each spectrum the nine parameters listed in Table [T] are extracted. These 
parameters may have contributions from different heights and/or solar features, and their 
response to solar activity is not identical. As described in Bertello et al. (2011), K3 
intensity, 0.5- A and 1-A emission indices are highly correlated with each other. Wavelength 
separation of V & R emission maximum does not correlate well with any other parameters, 
while K3 wavelength shows moderate correlation with most of the other parameters. 



3. Data Processing 

From December 2, 2006 to January 27, 2012 the ISS instrument has taken 2078 
measurements in the Ca II K spectral band. Each individual spectrum is normalized 
according to the procedure described in the previous section and the set of nine parameters 
listed in Table [1] is computed. Time series of these parameters are updated as new 
daily observations are taken and can be accessed from the SOLIS-ISS web site at 
http://solis.nso.edu/iss 
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3.1. Data Rejection 

Figure [1] sliows tlie temporal variation of some of the parameters investigated in 
this study. The other parameters exhibit a similar behavior, with the exception of Aks 
that does not show any significant temporal trend and was therefore excluded from our 
analysis. The rising phase of Cycle 24 is clearly visible in the temporal behavior of these 
parameters. The error bars of the individual measurements are typically smaller than the 
size of the points shown in the figure, and are omitted for clarity. Although the error bars 
for individual measurements are small, there is a small number of measurements which 
deviate significantly from the general population. These outliers were traced to specific 
(bad) observations. The first step in the data reduction consists of rejecting the outliers 
from the time series. This was achieved by fitting a cosine function with period of ~11 
years to the data and eliminating points that are 3-0" away from the model. The procedure 
was repeated until no points were eliminated. The number of points rejected depends on 
the particular time series, varying from 14 (Ik2V-Ik3)/(Ik2R-Ik3) to 44 (Ik2v/Ik3) out of 
2078 total points. 

3.2. Resampling and Filtering 

Each clear day the ISS takes typically 1-2 observations in the Ca II K bands, producing 
time series of parameters that are not evenly sampled in time. The problem of estimating 
the power spectral density from irregularly sampled data is in general more complicated 
than that of equidistant data. The direct Fourier transform and the method of Lomb-Scargle 
converge to the usual periodogram if the observations are equidistant, but suffer from a 
significant bias in the estimated spectra of unevenly sampled time series (Benedict et al. 
2000). Furthermore, the Lomb-Scargle periodogram (Scargle, 1982) is particularly designed 
to detect sinusoidal signals in white-noise unevenly time series. However, observational data 
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such as those investigated here often contain fractions of non-Gaussian noise or may consist 
of periodic signals with non-sinusoidal shapes. These properties of the data make the 
interpretation of the Lomb-Scargle periodogram more difficult and can lead to misleading 
estimates. Depending on the properties of the time series, many resampling methods are 
available (e.g., de Waele & Broersen 2000) that allow the use of more effective spectral 
estimators. Nearest neighbor resampling with the slotting principle (Broersen, 2009) is 
a technique that replaces an irregularly spaced sample by an equidistant signal with the 
resampling distance A. At every resampling node, the closest irregular observation is 
substituted if it is within half the slot width A from the resampling time. If no irregular 
observation falls within the slot width, the grid node is left empty as a missing observation. 
If the duty cycle of the resampled time series is reasonably high, missing observations 
can be filled via linear interpolation without introducing spurious peaks in the spectral 
estimator (de Waele & Broersen 2000). Because of the high temporal coverage of the ISS 
observations, a suitable value for the slot width is A = 1 day, with the resampling time 
set at Noon local time. This particular choice for the resampling time selects the daily ISS 
observation with the best signal-to-noise ratio. 

To eliminate spectral leakage from the very low frequency range, each time series was 
detrended using a 4-th order B-spline. Next, the signals were digitally filtered by means 
of a low-pass finite impulse response linear phase filter designed according to the recipe 
described in McClellan et al. (1979). The application of the filter before the spectral 
estimate is necessary to ensure that signal distortions due to aliasing are neghgible. The 
filter was designed to have a passband between and 0.2 cycles/day (5 days), a transition 
band between 0.2 and 0.4 cycles/day (2.5 days), and a stopband above 0.4 cycles/day. 
This was achieved using 21 coefficients, which gives an attenuation of ~65 db for the 
stopband and a maximum deviation of ~0.001 db in the passband. Finally, the data were 
pre- whitened to reduce the non-stationary components of the signal. 
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3.3. Spectral Estimation 

Among the available spectral estimators we selected the maximum entropy spectral 
estimator (MEM), which is equivalent to the autoregressive (AR) spectral estimator 
(Brockwell & Davis 1991). The maximum entropy spectral estimator is well suited for the 
analysis of our data because of its high spectral resolution, even with relatively short time 
series. If {xi} is a zero-mean Gaussian stochastic process A^(0,(T^), which models our time 
series, the spectral density S^u) at frequency u is given by 



|1 + Efc=i afcexp(-zA;z/)|2' 
where and p are the coefficients and order of the AR process: 

p 

Xi = '^akXi_k + Zi. (2) 

k=l 

Different estimation procedures exist as to how the AR coefficients, a^, and the white 
noise variance, Zi, are estimated in Equation [2] for any order p (e.g., Robinson & Treitel 
2000, Brockwell & Davis 1991, Percival & Walden 1993). The optimal choice of p needs to 
be considered carefully: if the number of poles, p, is too high, the method can introduce 
spurious and/or split peaks when applied to a noisy time series with a very large number 
of data points. A number of criteria have been suggested for selecting p (e.g., Brockwell & 
Davis 1991), but in most applications p is chosen empirically (Papoulis 1991). For our study 
we verified that as long as p is inside a reasonable interval, p G [A^/4, A^/2] (A^ being the 
number of data points), the same spectral features are found to be statistically significant 
for all values of p. A value of p = was used in our analysis. 

Unlike more traditional approaches based on Fourier techniques, the MEM lacks an 
easy and self-consistent procedure for evaluating the statistical significance of the spectral 
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estimates. A significance can be estimated with a permutation test (Good, 2000), which 
has two main advantages: 1) it is distribution-free, and 2) it is apphcable with any spectral 
estimator. The main disadvantage is that this method can be very computer-intensive, 
particularly with long time series. The implementation of the permutation test for spectral 
analysis requires the following steps: 

1. From the original sequence of data {x{n),n = 1,2, ... , N}, the power spectrum, S'(z/j), 
is computed for a set of frequencies {i^i}. 

2. A high number, M, of random sequences are obtained by the random permutation 
of the original data sequence, {x{n)}. For this work we used M = 50,000. Each 
sequence, {xj{n)}, obtained by random permutation has exactly the same data as the 
original sequence, but in a different order. 

3. M different power spectra, Sjiyi), are calculated from the M permutation sequences. 
The power spectrum is estimated using the same order p and for the same set of 
frequencies, {vi}, as the original power spectrum. 

4. The statistical significance of the spectral features of the original sequence is 
assessed for each frequency, Vi, by the achieved significance level (ASL) defined as 
ASL(z/j) = I{S{vi) > Sj{vi)}/M, where I{A} is the indicator function of event A 
(Pardo-Igiizquize & Rodriguez- Tovar 2005). I {A} is equal to 1 if event A happens 
and otherwise. Spectral values S{vi) with ASL(z/j) < a are considered statistically 
significant with confidence 100(1 - a)%, where the significance level a is chosen to be 
0.001 (99.9% confidence level) in this work. 

For comparison with the spectral estimation method described in this section, we 
computed the power spectrum of the full-length 1-A EM time series using the direct Fourier 
transform and the Lomb-Scargle periodogram. This comparison is shown in Figure [2l The 
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same re-sampled and interpolated time series was used for the Fourier and MEM analyzes, 
while for the Lomb-Scargle periodogram we used the one-a-day measurements (closest to 
Noon) and preserved the original time stamp. The frequencies at which the power was 
computed are the same for all three estimations. All three methods detect similar spectral 
features, but the MEM is able to isolate the ~ 27-day rotational component from the 
time series with a much higher signal-to-noise ratio. For example, within the investigated 
spectral range, the ratio between the primary and secondary peaks is about a factor five 
for the MEM while it is significantly lower for the other two estimators. In addition, the 
Lomb-Scargle periodogram clearly shows a much more complex spectral structure compared 
to the other two estimates. 



4. Results and Discussion 

The ISS parameter time series mostly cover the extended minimum of Cycle 23/24 
when the activity was low, making the detection of SDR more difficult. The main purpose of 
this study is to show that, under these less than favorable conditions, the selected technique 
of spectral analysis plays a substantial role in the ability to detect the signature of SDR. 
The major constraint imposed by the low level of solar activity covered by this investigation 
is in the length of the temporal window that needs to be selected for the spectral estimator. 
As discussed below, under these circumstances, at least 900 days of ISS data are required 
to provide a (statistically) significant detection of SDR. As the ISS program continues, the 
observations will cover the full duration of Cycle 24, and it is expected that shorter time 
series will be needed. 
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4.1. Spectral Analysis 

We present here results from the spectral analysis performed on the full-length 
time series and on the first and last 2/3 portions of the ISS data. To explore the time 
evolution of individual components in the power spectrum we also subdivided every time 
series in overlapping segments, shifted by 1 day, and computed the spectra of these 
segments (Section 4.2). In addition, we applied the same analysis to a time series of 
disk-averaged longitudinal magnetic field flux density (MF) measurements derived from 
daily SOLIS Vector Spectromagnetograph (VSM) magnetograms taken in the Fel 630.15 
nm spectral line (Pietarila et al. 2012). The addition of an independent time series is 
necessary to eliminate the possibility that our choice for the spectral estimator is biased 
towards a particular data set. The selection of the MF time series was made because it 
provides the same time coverage as the ISS observations, and also because the emission 
parameters derived from the Ca K spectral line are closely related to the unsigned magnetic 
flux (e.g., Ortiz Sz Rast 2005). The MF data are available from the SOLIS website at 



http://solis.nso.edu/vsm/vsm_mnfield.html It is important to emphasize here that the 
power spectra of individual time series are not expected to show the same peaks at the 
same level of significance. There are two main reasons for that: 1) different parameters are 
computed from different portions of the Ca 11 K line profile, and therefore have different 
sensitivity to solar activity; 2) parameters computed from line-intensity values are more 
affected by possible calibration issues than those computed from wavelength separations 
(see Tabled] for a description of the parameter used in this work). 

Figure [3] shows the power spectra of the full-length MF and four of the eight ISS 
full-length parameter time series investigated in this study. All time series were divided by 
their standard deviation to get a zero mean unit variance time series. The power spectrum 
of the MF time series is shown on the top, followed by the two spectra corresponding to 



- 14 - 



parameters computed from the wavelength separations of Ca K hne features. The bottom 
two spectra are Ca K intensity-related parameters. This Figure clearly shows the different 
features associated to the two groups of parameters, and it is also representative for the 
other parameters listed in Table [H As shown in column 2 of Table [1], in eight of the nine 
time series investigated in this work the tallest peak above the 99.9% confidence level 
corresponds to a periodicity between 27.4 and 28.0 days. These numbers are within the 
range of values expected for the solar rotation period. In general, parameters defined in 
terms of wavelength separation typically show a spectrum with a very predominant peak, 
while intensity-related parameters show spectra with a more complex structure. The only 
parameter showing a somewhat more erratic power spectrum is Ak2R — Ak2V, the wavelength 
separation of the two emission maxima. The power spectrum of this time series has a 
predominant peak corresponding to a rotational period of 141 days, most likely the result 
of low sensitivity of the parameter to solar activity. 

The detection of solar rotational periods in the power spectra of Ca II K parameter 
time series can be affected by unrelated signatures of other processes: sunspot evolution 
(typically weeks), faculae evolution (up to ~60 days), and active longitudes characterized 
by the frequent emergence of new magnetic flux (up to ~200 days). These processes can 
reduce the coherence of the rotational signal to the point where it is no longer detectable, 
and some time series can be affected more than others. We believe this to be the case for 
the power spectrum of Ak2r — Ak2v- In fact, a recent analysis (Pevtsov & Bertello 2012) 
of basal and active sun emission components of the CaK line profiles indicates that the 
wavelength separation between the active Sun components does not change with the phase 
of the solar cycle. 

We test the ability of our method to detect SDR during the solar cycle by computing 
the power spectra of the first and last 2/3 portions of the individual time series. The 
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results are shown in Figure ID In the left column are the spectra computed from the first 
2/3 portion of the data, while in the right column we show the spectra from the last 2/3 
portion of the time series. The second portion of the data produces spectra characterized 
by the presence of a single predominant peak, while the spectra of the first 3.5 years of data 
show a more complex structure. This is not surprising since the first 2/3 of the data mainly 
cover the extended minimum of Cycle 23, and therefore the signature of rotation is not as 
dominant. A comparison between the two sets of spectra shows also a clear difference in the 
location of the tallest peak: around ~26.3 days for the first 3.5 years of data, consistent with 
activity near the equator (~ 7°, using the formula in Snodgrass & Ulrich 1990), and ~27.7 
days for the portion of the time series after August 2008. This distinct period of ~27.7 days 
is in agreement with features rotating at higher latitude bands (~ 30°), as expected from 
the transition between Cycles 23 and 24. The power spectrum of the first portion of the 
MF 630.15 nm time series shows two very predominant peaks, located at 27.4 days and 26.3 
days. While the secondary peak at 26.3 days is in very good agreement with the results 
from the ISS time series for the same time frame, the tallest peak at 27.4 days is more 
consistent with the ISS results from the second part of the data. A possible interpretation is 
that late in the time series covered by the first segment of the data, a new-cycle of activity 
was already present at higher latitudes. These high latitude fields were relatively weak, 
and did not contribute significantly to the disk-integrated Ca K profiles. Figure [5] shows 
the relationship between the VSM total unsigned flux and the 1-A EM index for three 
different periods. The correlation during the raising phase of Cycle 24 (2010-2012) is strong 
and statistically significant, while during the dechning phase of Cycle 23 (2007-2008) the 
correlation is weak, although still statistically significant. During solar minimum (2009) no 
correlation is found. Therefore, the Ca K emission is a good proxy for magnetic activity 
only for periods of intermediate and high levels of solar activity, which indirectly supports 
our above interpretation that the presence of weak magnetic fields at high latitudes may 
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explain the difference in periodicity between the Ca K and MF results. 

4.2. Time-frequency Analysis 

To further confirm the results shown in Figure HI we have performed a time-frequency 
analysis to investigate the time evolution of the individual components of the power 
spectra computed from the different time series. We have subdivided every time series in 
overlapping segments of 900 days, with a 1-day time shift between consecutive segments, 
and computed the power spectra according to the recipe described in the previous section. 
According to our tests with temporal windows of different sizes, a 900-day sliding window 
provides the best compromise between time and frequency resolution. Figures [6HH] show 
time-frequency plots computed for two of the ISS parameter time series and for the MF 
data. These plots are representative of the other time series. The signature of SDR is clearly 
visible in these plots: most of the contribution to the power spectral density of the signal is 
concentrated into a narrow frequency band, whose central value varies in time in the range 
0.035 to 0.038 cycles/day (period range from 28.6 to 26.3 days). This change in period over 
time is consistent with the migration of activity from latitudes bands at ~ 35 degrees to 
the equator. To large degree, this curve is remarkably continuous and smooth at all times. 
These figures also show, for the ISS measurements, some additional power distribution over 
different time-frequency regions. This is likely the result of some of the physical processes, 
e.g., sunspot and facular evolution, discussed in the next section, and possible artifacts 
related to data calibration. A close look at the time-frequency distributions around the 
period of solar rotation shows some differences in the temporal behavior of the power 
spectra. For example, the rotational period during the declining phase of Cycle 23 deduced 
from the MF data seems to be slightly longer, and more persistent, than the one derived 
from the ISS measurements. This could be explained by changes in the correlation between 
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magnetic field and Ca K proxies with activity level (Figure |5]). 

4.3. Uncertainty in the Rotational Period 

The detection of SDR from disk-integrated measurements is affected not only by the 
physical processes mentioned above but also by the particular analysis performed on the 
data. With the approach described in this paper, the main factors limiting the identification 
and accuracy of the rotational components in the power spectrum are: (1) the resampling 
and gap filling of the signal, (2) the choice of the order p of the AR process used to model 
the time series, and (3) the length of the time series window. The effect of the resampling 
and gap filling has been tested by computing the power spectrum of the window function. 
We did not find any significant contribution from the window function to the power spectral 
density of the signals in the region of interest, between ~ 22 and ~ 33 days. The order 
p mainly changes the width of the rotational peaks (higher values of p produce narrower 
rotational peaks) without changing (significantly) their central frequency. We can determine 
the uncertainty in the deduced rotational period from the width of the peaks in the power 
spectrum and from the variation in their central location as a result of small changes (10%) 
in the selected order p of the AR process. Using this approach, the estimated uncertainty 
of the values hsted in Table 1 is approximately ±0.3 days. 

4.4. Length of the Time Series 

The choice for the length of the time series window used in the analysis described in 
section 3 depends on three major factors: 1) the nature and properties of the data, 2) the 
level of activity, and 3) the period of rotation. For the case of the Sun, a slowly rotating 
star, our results using ISS data show that a 900-day window was necessary to provide a 
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statistically significant detection of SDR over the whole period covered by the observations, 
which coincide with a prolong period of extremely low activity. Shorter time series (few 
hundred days) may affect the statistical significance of the peaks and impact the overall 
spectral distribution by producing spectra with a distribution of peaks over a broad spectral 
range. Furthermore, the peaks may often be inconsistent from one observing window to 
another, making it very difficult to disentangle contributions from different sources. Using 
chromospheric full-disk Ca II K SPO data in 200-day intervals, Donahue & Keil (1995) 
were only able to detect rotation in 42% of the observing windows, mostly during periods 
of high solar activity. 

To investigate the effect of the activity level on the detectability of solar rotation on 
shorter time intervals, we tested both the Lomb-Scargle periodogram and the MEM on a 
sequence of 250-day ISS data segments. The MEM performed significantly better than the 
Lomb-Scargle periodogram, producing spectra of the same quality as those shown on the 
right panel of Figure S] with a clear detection of the rotational peak for about 38% of the 
observing windows during periods of both medium and high solar activity. The determined 
rotation periods from these spectra are consistent with the results obtained using a 900-day 
sliding time window. This test shows that the MEM can effectively be used with relatively 
short time series, provided there is enough activity for the signature of rotation to be 
present in the signal. 

The application of this method to stellar measurements is more complex, but we see no 
reasons this technique would not work effectively also for the stellar case. Quasi-continuous 
Ca II K and H stellar measurements from ground-based observatories are typically limited 
in duration to a maximum of about 6 months. Our choice of 900-day and 250-day windows 
was mainly constrained by the fact that the Sun is a slowly rotating star and the ISS 
observations covered the extended solar minimum, with periods of no plages or other weak 
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remnants of active regions on the Sun. 

Young solar-like stars, however, are fast rotators with periods of a few days. Thus, 
for example, for a star with 5 days rotation period and low spot-activity it should require 
about (900/27) -5 = 167 days of observing time with a 58% duty cycle. For a star with high 
spot-activity this number would be (250/27) -5 = 46 days. The above estimates refer to a 
detection of the rotation with a high significant level of confidence; the rotation modulation 
can be detected with even shorter observing windows albeit at a lower confidence level. This 
scenario is slightly different for the case of stellar light curve time series. The rotational 
modulation in light curves is caused by dark spots as well as bright faculae. In slowly 
rotating stars with relatively low level of activity, such as the Sun, these two contributions 
merge together to reduce the sensitivity of the signal to rotation modulation. On the other 
hand, dark spots dominate active and young stars while faculae are predominant in weakly 
active stars during activity minimum. In these two cases the signature of rotation can be 
observed from light curves if a sufficiently high signal-to-noise ratio is achieved. Several 
ground- and spaced-based projects provide light curves over months or years that can be 
investigated using the proposed technique. In particular, with a lifetime of at least 3.5 years 
and nearly continuous observations, the Kepler mission (Borucki et al. 2010) is delivering 
stellar light curves of unprecedented time sampling and photometric precision. 

4.5. Simplified Flux Transport Model 

To further validate our selection of parameters for spectral analysis, and to investigate 
the effects of properties of active Sun features on the detectability of SDR, we synthesized 
artificial data using a simplified flux transport model. The model evolves the radial 
magnetic field by the effects of flux emergence, differential rotation, meridional flow, and 
diffusion. First, we created a magnetogram representing the entire solar surface covering 
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360 degrees in longitude and 180 degrees in latitude. The magnetogram is made of two 
components representing strong and weak magnetic fields. For weak fields we use a random 
distribution of small (in flux and size) bipolar elements. To simulate stronger fields of 
active regions, we insert larger bipolar pairs in selected latitudes and longitudes. The 
magnetogram is evolved by applying the known sunspot profile of SDR (Newton & Nunn 
1951), diffusion, and meridional circulation. Throughout the simulations, the magnetogram 
is updated by inserting new active regions. The location (latitude and longitude) of 
simulated active regions, their size and rate of emergence are taken from actual sunspot 
observations taken at the Kislovodsk Solar Station of Pulkovo Observatory in Russia 
(http://www.solarstation.ru/?lang=en&id=lastdata) during the same period as the ISS 
observations analyzed in this study (December 2, 2006- January 27, 2012). The magnetic 
flux of simulated active regions is defined using the known relation (Houtgast & van Sluiters 
1948) between sunspot area (S, in millionth of solar hemisphere) and magnetic flux ($ in 
Gauss): 

^ 3700 X S 
$ = . 

S + 66 

Since the Kislovodsk dataset does not contain information about the tilt or polarity 
orientation of active regions, for the modeling we assume that all low latitude (< 20 degrees) 
active regions observed prior to year 2009 (approximate time of solar minimum) are oriented 
in agreement with the Hale polarity rule for solar Cycle 23. High latitude active regions 
observed prior to 2009 and all regions observed from 2009-2012 are assigned a polarity 
orientation in correspondence with the solar Cycle 24. For simplicity, all active regions are 
assigned the same tilt relative to the equator. We experimented with tilts following Joys 
law (Hale et al. 1919) but found no significant differences in the results. New active regions 
are "inserted" in the simulated magnetogram based on the time of the first observation 
of the region. As an additional simplification, emergence of multiple active regions on the 
same day was ignored (only one active region was inserted in this limited number of cases). 
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Unlike in more realistic flux transport models (e.g., Jiang et al. 2011), diffusion in our 
model is done by convolving the magnetogram with an inverse Gaussian filter. 

To investigate the role of diffusion in deriving solar rotation from Sun-as-a-star data, 
we run the model with six spatial filters corresponding to different parametrizations for 
diffusion (width d of Gaussian filter): d=5, 10, 15, 30, 40, and 50. We found that d=30 
results in the complete disappearance of an active region with strong magnetic flux over 
approximately 5- 6 solar rotations, thus providing a proxy for a reasonable solar diffusion 
rate. Based on this, we refer to models run with d=5, 10, 15 of "fast" diffusion, and 

d=40, 50 clS CclSGS of "slow" diffusion. As a qualitative reference, with d=5 newly injected 
active region magnetic field completely disappears in 2-3 days, and with d=50 the presence 
of diffusing flux is seen for 10+ solar rotations. 

Using the set of simulated 360 degree x 180 degree maps, we create full disk line-of-sight 
magnetograms for each day of observations during the ISS observing period. These daily 
magnetograms are used to create a time series of Sun-as-a-star unsigned and signed fluxes. 
Each series covers 1883 days (starting on December 2, 2006 and ending on February 27, 
2012) and has a uniform time cadence of one day. The simulated time sequences are 
re-sampled as the ISS observations and subjected to the same data analysis. 

Examples of time-frequency plots computed for the cases of "slow" (d = 40) and 
"fast" (d = 5) diffusion are shown in Figures [9] and [TOl The case of "fast" diffusion has 
no clear signature of solar rotation, indicating that the fast diffusing features do not live 
long enough to be detected by our analysis. This case of "fast" diffusion may explain 
previous failed attempts in detecting differential rotation from disk integrated white-light 
spectra. Although the white-light observations have a contribution from features (such as 
sunspots) that follow solar rotation, these features decay rapidly (to small pore sizes) and 
their contribution to the white-light disk-integrated signal is diminished. 
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When the diffusion time increases (Figure [TOl) . our technique is able to unambiguously 
detect the signature of differential rotation. The excellent agreement between the results 
from the observed and simulated data is a further validation of the method described in 
this paper. In addition, we also experimented with randomly injecting magnetic bipoles 
with a varying rate of emergence. We found that even when the active region emergence 
rate is low (up to one active region in 300+ days), we were able to detect the solar rotation 
for normal (d=30) and slow diffusion rates. For rapid diffusion (d=5), we found that the 
rotation rate is not well detectable even at very high rate of emergence of new active regions. 
To investigate the effect of active longitudes (when active regions emerge preferentially at 
two longitudes separated by about 180 degrees), we run our code with randomly emerging 
active regions near these two preferred longitudes. Our analysis shows no definite splitting 
in rotation frequencies, which would indicate the presence of active longitudes. 

5. Conclusions 

In this paper we show that even under unfavorable conditions of very low levels of 
solar activity it is possible to detect what can be characterized as the signature of surface 
differential rotation from disk-integrated measurements. We present evidence of strong 
components in the power spectra computed from time series of parameters derived from the 
temporal variations of Ca II K ISS line profiles and the MF time series that are consistent 
with the rotation rate of active regions migrating on the solar disk during the declining 
phase of Cycle 23 and the rising phase of Cycle 24. The length of the time series window 
used in our analysis is significantly longer than the one adopted in previous studies, but 
essential to characterize the pattern of differential rotation during periods of low magnetic 
activity. Longer time series also have the advantage of producing power spectra with a 
reduced number of significant peaks over a large spectral band, making it easier to identify 
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the different contributions to the signal. In our tests, the adopted spectral estimator has 
been proven to be significantly more effective than the Lomb-Scargle periodogram method 
in providing spectra with a reduced number of significant peaks. This is very important 
when prior knowledge of the rotation rate profile is uncertain, as is the case for stellar 
observations. As some methods of spectral analysis may produce multiple peaks in a broad 
range of periods, it is important to correctly identify the peak corresponding to a true stellar 
rotation. The results from the ISS observations and from our numerical model suggest that, 
in addition to the significance level, the continuity of a specific frequency in power spectra 
can be used as a criterion. For example. Figure M shows the presence of several peaks of 
different duration, but the peak around 27-28 days is the only one that persists throughout 
the entire data set. 

The results of our numerical modeling suggest that the diffusion rate of active regions 
plays the most important role in detection of solar rotation from Sun-as-a-star observations. 
Although we use a simplified flux transport model, our results show the significance of the 
diffusion time-scale parameter of active regions in the context of Sun-as-a-star differential 
rotation, with possible ramifications to stellar applications as well. The rate of emergence 
and the presence of active longitudes seem to play a less relevant role. Due to the 
simplifications in the adopted transport model, our results are mostly qualitative in nature. 
An important challenge therefore is to quantitatively assess the implications of diffusion 
on determining SDR from disk-integrated measurements using more realistic flux transport 
models, a study that is beyond the scope of this paper. 

Although some differences exist in the results obtained from the various time series 
investigated in this study, the overall excellent agreement in the detection of SDR is a 
strong validation of the adopted method. The use of simulated data, with a known rotation 
profile, is an important confirmation of the quality of our methodology. We believe that this 
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approach can be successfully applied to other time series, particularly stellar measurements, 
even when the duty-cycle is not very high. In fact, our re-sampled ISS time series have a 
duty cycle of approximately 58%. 

SOLIS data used here are produced cooperatively by NSF/NSO and NASA/LWS. The 
NSO and NOAO are operated by the Association of University for Research in Astronomy, 
Inc. (AURA), under cooperative agreements with the National Science Foundation. 
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Table 1: Measured synodic rotation periods calculated from the power spectra of the eight ISS 
Ca II K parameter time series, and SOLIS-VSM Fel 630.15 nm mean longitudinal magnetic 
field flux density (MF). Rotational periods have been calculated from spectra of the full- 
length time series (column 2), and from spectra computed using the first (column 3) and the 
last (column 4) 2/3 portions of the time series. Only the period of the tallest peaks above 
the 99.9% confidence level is hsted. The estimated uncertainty of these values is ±0.3 days. 



Parameter 


Full-length 

Rotat: 


First 2/3 

ion ])eri()(l (d 


Last 2/3 


1-A EM 


27.7 


26.3 


27.8 


0.5-A EM 


27.7 


26.3 


27.7 


Ik3 


27.8 


26.3 


27.8 


Wilson-Bappu 


27.4 


26.4 


27.6 


Akir — Akiv 


27.7 


26.2 


27.7 


(Ik2V-Ik3)/ (Ik2R-Ik3) 


27.5 


26.0 


27.7 


Ak2R — Ak2V 


141.0 


none 


none 


Ik2v/Ik3 


27.7 


51.3 


27.8 


MF 630.15 nm 


28.0 


27.4 


27.8 



Note. — 1-A EM is the 1-A emission index, 0.5-A EM is the 0.5-A emission index, Iks is the intensity 

in the line core, Wilson-Bappu is the wavelength separation between outer edges of K2, Akir — Akiv is the 
wavelength separation of V & R emission minimum, (Ik2V-Ik3 ) / (Ik2R-Ik3 ) is the line asymmetry, Ak2R, ^ Ak2V 
is the wavelength separation of V & R emission maximum, and Ik2V /Ik3 is the ratio of the blue emission 
peak to the intensity in the line core. The violet and red wings of the line profile are indicated by V and R 
respectively. 
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Fig. 1. — Four of the eight ISS Ca II K parameter time series listed in Table 1. The 
period of time covered by these observations is from December 2, 2006 to January 27, 
2012. Daily updates of all parameter time series are available from the ISS web site at 



http: //solis. nso.edu/iss 
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Period (days) 

200.0 100.0 66.7 50.0 40.0 33:3 28.6 25.0 22.2 20.0 18.2 16.7 




0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 

Frequency (cycles/day) 

Fig. 2. — Comparison of three different spectral estimators applied to the full-length 1-A EM 
time series. The same re-sampled and interpolated data was used for both the Fourier and 
MEM analyzes, while for the Lomb-Scargle periodogram we used the one-a-day measure- 
ments (the closest to Noon) and preserved the original time stamp. The frequencies at 
which the power was computed are the same for all three estimations. Each power spectrum 
is normalized to its maximum value. 



- 31 - 



Period (days) 

200.0 100.0 66.7 50.0 40.0 33.3 28.6 25.0 22.2 20.0 18.2 16.7 



CO 

o 

CD 
Q. 
C/3 

T3 
CD 
_N 

"cO 

E 




0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 

Frequency (cycles/day) 

Fig. 3. — Spectral estimation of the four full-length time series shown in Figure [1] and 
the SOLIS-VSM Fel 630.15 nm mean longitudinal magnetic field flux density (MF). The 
horizontal dashed-line in each spectrum indicates the 99.9% confidence level. The list of 
rotation periods corresponding to the tallest peak above the 99.9% confidence level is given 
in Table [H which also contains the results from the other time series investigated in this 
study. 
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Period (days) Period ( , ^ 

100.0 66.7 50.0 40.0 33.3 28.6 25.0 22.2 20.0 18.2 16.7 100.0 66.7 50.0 40.0 33.3 28.6 25.0 22.2 20.0 18.2 16.7 




1 - 




1 - 



8/2008- 1/2012- 



0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 
Frequency (cycles/day) Frequency (cycles/day) 



Fig. 4. — Same as in Figure [31 but the power spectra were computed using the first (left) 
and the last 2/3 (right) portions of the full-length time series. The time interval covered 
by each segment is indicated in the plots. The horizontal dashed-line indicates the 99.9% 
confidence level. The list of rotation periods corresponding to the tallest peak above the 
99.9% confidence level is given in Table [H which also contains the results from the other 
time series investigated in this study. 
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Fig. 5. — Comparison of VSM total unsigned flux and ISS 1-A EM index for the the 
dechning phase of Cycle 23 (left), solar minimum (center), and raising phase of Cycle 24 
(right). Shown in sohd line is the computed linear fit. The Pearson correlation coefficient, r, 
is given in the upper left corner of each panel. We computed the t-value for the coefficient, 
r, to assess the statistical significance of each correlation. The results are (for the left to 
the right panels): 13.6 (2.59), 2.48 (2.69), and 37.3 (2.59). Indicated in parenthesis are the 
expected 99% cutoff values in the corresponding two-tailed Student's t distribution. 
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Fig. 6. — ISS Ca II K3 core intensity time-frequency distribution from observations taken 
during the period 12/2/2006-1/27/2012. A 900-day sliding time window was used, with a 
difference between consecutive segments of 1 day. The beginning of each segment, in days 
since 2/12/2006, is indicated on the left y-axis with the year shown on the right y-axis. The 
power spectral density above/below the 99.9% confidence level (CL) is indicated by the color 
bar. 
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Fig. 7. — Same as Figure [6] but for the ISS Akir — Akiv parameter. 
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Fig. 8. — Same as Figure O but for the SOLIS-VSM Fel 630.15 nm mean longitudinal 
magnetic field flux density time series (MF). 
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Fig. 9. — Same as Figure |6l but for artificial data from a simplified flux transport model 
with "fast" diffusion (d=5, see text). 
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Fig. 10. — Same as Figure |9l but for the case of "slow" diffusion (d=40). 



